﻿rm(list = ls())
library(tidyverse)

data.table::fread("/Pub/Users/wangyk/project/Poroject/0524_RS10_HNSC_ARGININE_AND_PROLINE_METABOLISM/output/forestplot_univ_multi_cox_GSE65858_OS/multiCox_forestplot_GSE65858_OS_res.txt") %>% 
as_tibble()

load("/Pub/Users/wangyk/project/project_test_Pre_Sales/P220112001_SKCM_3+/res/major_res.RData",verbose = T)

major_res %>% names()

pheno_1 <- major_res$training %>% mutate(Stage = factor(Stage),Age = factor(Age))

source("/Pub/Users/wangyk/Project_wangyk/Codelib_YK/some_scr/univ_multi_Cox_forestplot_clinical_feature.R")
load("/Pub/Users/wangyk/project/Poroject/F210823005_非小细胞肺癌血小板/out/forest_data.RData")
head(forest_data)
colnames(forest_data) 

c("Smoking","EGFR","ALK_eml4","KRAS","Histological.Type","Score")


a <- univ_multi_Cox_forestplot(
    clinical = forest_data, saveplot = F, out_dir = "./",
    features_col = c("Smoking", "EGFR", "ALK_eml4", "KRAS", "Histological.Type", "Score"),
    var_name = NULL, width = 7.5, points_color = "#005670e3", multi_after_univcox = F
)

library(gtsummary)
library(survival)

# build survival model table
t2 <-
  coxph(Surv(ttdeath, death) ~ trt + grade + age, trial) %>%
  tbl_regression(exponentiate = TRUE)

Forest_DF <- function(clinical = NULL, features_col = c("Smoking", "EGFR", "ALK_eml4", "KRAS", "Histological.Type", "Score"),
                      SaveFile = T, od = "./", w = 5, h = 3,p_cut = 2,var_name = NULL) {

    require(tidyverse)
    require(survival)
    require(survminer)
    require(cowplot)
    require(ggpubr)
    require(flextable)
    require(officer)
    require(gtsummary)

    if(is.null(var_name)){
        var_name <- Sys.time() %>% as.character()
    }

    levels <- map(features_col, \(x){
        clinical %>%
            select(any_of(x)) %>%
            pull() %>%
            factor() %>%
            as.list() %>%
            lvls_union()
    })
    names(levels) <- features_col

    n_df <- map(seq_along(features_col), \(x){
        clinical %>%
            select(features_col[x]) %>%
            table() %>%
            as_tibble() %>%
            rename("levels_name" = 1)
    })
    names(n_df) <- features_col

    features_chr <- map(features_col, \(x){
        str_c(x, " ", paste0("(", levels[[x]][-1], " vs. ", levels[[x]][1], ")"))
    }) %>%
        unlist() %>%
        as_tibble()

    n_chr <- map(features_col, \(x){
        n_df[[x]] <- n_df[[x]] %>% mutate_if(is.numeric, as.character)

        n_ref <- n_df[[x]] %>%
            filter(levels_name == levels[[x]][1]) %>%
            pull(n)

        n_levels <- map_chr(levels[[x]][-1], function(y) {
            n_df[[x]] %>%
                filter(levels_name == y) %>%
                pull(n)
        })

        str_c(paste0(n_levels, " vs. ", n_ref))
    }) %>%
        unlist() %>%
        as_tibble()

    vars_for_table <- clinical %>%
        dplyr::select(any_of(features_col)) %>%
        colnames(.)
    nv_var <- vars_for_table

    univ_formulas <- sapply(vars_for_table, function(x) as.formula(paste("Surv(time, status)~", x)))
    univ_models <- lapply(univ_formulas, function(x) coxph(x, data = clinical))

    univcox_clinical_res <- map_dfr(.x = univ_models, .f = function(x) {
        result <- summary(x)

        result2 <- data.frame(
            var = rownames(result$coefficients),
            coef = result$coefficients[, 1],
            Hazard_Ratio = result$coefficients[, 2],
            se = result$coefficients[, 3],
            # p.value = result$coefficients[, 5],
            p.value = result$sctest["pvalue"],
            lower_.95 = result$conf.int[, "lower .95"],
            upper_.95 = result$conf.int[, "upper .95"],
            logrank_pvalue = result$sctest["pvalue"],
            wald_pvalue = result$waldtest["pvalue"],
            Likelihood_pvalue = result$logtest["pvalue"]
        )

        result2 <- result2 %>%
            mutate(HR = paste0(
                sprintf("%.2f", Hazard_Ratio),
                "(",
                sprintf("%.2f", lower_.95),
                "-",
                sprintf("%.2f", upper_.95),
                ")"
            ))
    })

    rownames(univcox_clinical_res) <- univcox_clinical_res$var
    univcox_clinical_res$var <- features_chr$value

    univ_cox_df <- univcox_clinical_res %>%
        mutate(N = n_chr$value) %>%
        mutate(P = ifelse(p.value < 0.001, "<0.001", sprintf("%.3f", p.value))) %>%
        select(var, N, HR, P) %>%
        rename("HR(95%CI)" = HR)

    p_univ_cox <- univ_cox_df %>%
        ggtexttable(., rows = NULL, theme = ttheme("blank")) %>%
        tab_add_hline(at.row = 1:2, row.side = "top", linewidth = 2) %>%
        tab_add_hline(at.row = dim(univ_cox_df)[1] + 1, row.side = "bottom", linewidth = 2) %>%
        tab_add_title(
            text = "Univariate  Cox",
            face = "bold",
            size = 16,
            padding = unit(1.5, "line"),
            just = "left"
        )

    # MultiCox Part-------------

    # p_cut <- .05
    do_feature_select <- function(p_cut) {
        col_names_left <- map(univ_cox_df %>% filter(P < p_cut) %>% pull(1) %>% str_split(" \\("), ~ .x[[1]]) %>%
            unlist() %>%
            unique()

        vars_for_table <- clinical %>%
            dplyr::select(any_of(col_names_left)) %>%
            colnames(.)

        levels <- map(col_names_left, \(x){
            clinical %>%
                select(x) %>%
                pull() %>%
                factor() %>%
                as.list() %>%
                lvls_union()
        })
        names(levels) <- col_names_left

        features_chr <- map(col_names_left, \(x){
            # x <- 'TCGA_Subtype'
            str_c(x, " ", paste0("(", levels[[x]][-1], " vs. ", levels[[x]][1], ")"))
        }) %>%
            unlist() %>%
            as_tibble()

        return(list("features_chr" = features_chr, "col_names_left" = col_names_left))
    }

    features_new <- do_feature_select(p_cut)

    if (nrow(features_new$features_chr) < 2) {
        features_chr <- features_chr
        vars_for_table <- vars_for_table

        warning(str_glue("在p小于 {crayon::blue(p_cut)} 的时候，只有一个临床信息保留，所以使用全部变量进行多因素Cox分析。"))
    } else {
        features_chr <- features_new$features_chr
        vars_for_table <- features_new$col_names_left
    }

    clinical_multi <- clinical %>%
        dplyr::select(any_of(vars_for_table)) %>%
        na.omit()

    multicox_formulas <-
        as.formula(paste(
            "Surv(time, status)~",
            paste0(sep = "`", vars_for_table, sep = "`", collapse = "+")
        ))
    # top 4 colunm：sample,time,status,riskgroup

    multicox_cox_res <- coxph(
        formula = multicox_formulas,
        data = clinical
    )

    multicox_cox_res_summary <- summary(multicox_cox_res)

    multicox_cox_res_df <- data.frame(
        var = rownames(multicox_cox_res_summary$coefficients),
        coef = multicox_cox_res_summary$coefficients[, 1],
        Hazard_Ratio = multicox_cox_res_summary$coefficients[, 2],
        se = multicox_cox_res_summary$coefficients[, 3],
        p.value = multicox_cox_res_summary$coefficients[, 5],
        lower_.95 = multicox_cox_res_summary$conf.int[, "lower .95"],
        upper_.95 = multicox_cox_res_summary$conf.int[, "upper .95"],
        logrank_pvalue = multicox_cox_res_summary$sctest["pvalue"],
        wald_pvalue = multicox_cox_res_summary$waldtest["pvalue"],
        Likelihood_pvalue = multicox_cox_res_summary$logtest["pvalue"]
    )

    multicox_cox_res_df <- multicox_cox_res_df %>%
        mutate(HR = paste0(
            sprintf("%.2f", Hazard_Ratio),
            "(",
            sprintf("%.2f", lower_.95),
            "-",
            sprintf("%.2f", upper_.95),
            ")"
        ))

    multicox_cox_res_df$var <- features_chr$value

    multi_cox_df <- multicox_cox_res_df %>%
        mutate(P = ifelse(p.value < 0.001, "<0.001", sprintf("%.3f", p.value))) %>%
        select(var, HR, P) %>%
        rename("HR(95%CI)" = HR)

    p_multi_cox <- multi_cox_df %>%
        ggtexttable(., rows = NULL, theme = ttheme("blank")) %>%
        tab_add_hline(at.row = 1:2, row.side = "top", linewidth = 2) %>%
        tab_add_hline(at.row = dim(multi_cox_df)[1] + 1, row.side = "bottom", linewidth = 2) %>%
        tab_add_title(
            text = "Multivariate Cox",
            face = "bold",
            size = 16,
            padding = unit(1.5, "line"),
            just = "left"
        )


    p_all <- ggpubr::ggarrange(p_univ_cox, p_multi_cox, align = "h")

    # title <- ggdraw() + draw_label("TCGA", fontface='bold' ,hjust = 0, vjust = 0) #nrow(univ_cox_df)
    # ggarrange(title,p_all,ncol = 1,height = c(.1,1))


    # new_style part-----------
    tbl_uv <- tbl_uvregression(
        clinical[c(nv_var, "time", "status")],
        method = coxph,
        y = Surv(time, status),
        exponentiate = T,
        pvalue_fun = function(x) ifelse(x < 0.001, "<0.001", sprintf("%.3f", x))
    )

    tbl_mv <- coxph(
        formula = multicox_formulas,
        data = clinical
    ) %>%
        tbl_regression(
            exponentiate = TRUE,intercept = T,
            pvalue_fun = function(x) ifelse(x < 0.001, "<0.001", sprintf("%.3f", x))
        )%>%
        add_n()

    tbl_merge_uv_mv <- tbl_merge(
        tbls = list(tbl_uv, tbl_mv),
        tab_spanner = c("**Univariate Cox**", "**Multivariate Cox**")
    )

    # 设置word格式
    # sect_properties <- prop_section(
    #     page_size = page_size(
    #         orient = "portrait",
    #         width = 21 / 2.54, height = 29.7 / 2.54
    #     ),
    #     type = "continuous",
    #     page_margins = page_mar()
    # )

    if (isTRUE(SaveFile)) {
        if(!dir.exists(od)){
            dir.create(od,recursive = T)
        }

        walk(c("p_univ_cox", "p_multi_cox"), function(x) {
            ggsave(plot = get0(x), filename = str_glue("{od}/Figure_{x}.pdf"), width = w, height = h)
        })

        ggsave(plot = p_all, filename = str_glue("{od}/Figure_Univ_Multi_Cox.pdf"), width = w * 1.75, height = h)

        walk(c("univ_cox_df", "univcox_clinical_res", "multicox_cox_res_df", "multi_cox_df"), function(x) {
            write.table(x = get0(x), file = str_glue("{od}/Table_{x}.txt"), quote = F, sep = "\t")
        })

        tbl_merge_uv_mv %>%
            as_hux_xlsx(file = str_glue("{od}/Table_Clinical_Univ_Multi_Cox_{x}.xlsx"))
        
        tbl_uv %>%
            as_hux_xlsx(file = str_glue("{od}/Table_Clinical_Univ_Cox_{x}.xlsx"))

        tbl_mv %>%
            as_hux_xlsx(file = str_glue("{od}/Table_Clinical_Multi_Cox_{x}.xlsx"))        


        tbl_merge_uv_mv %>%
            as_flex_table() %>%
            flextable::save_as_docx(
                path = str_glue("{od}/Table_Clinical_Univ_Multi_Cox_{x}.docx")
                # pr_section = sect_properties
            )
        tbl_uv %>%
            as_flex_table() %>%
            flextable::save_as_docx(
                path = str_glue("{od}/Table_Clinical_Univ_Cox_{x}.docx")
                # pr_section = sect_properties
            )
        tbl_mv %>%
            as_flex_table() %>%
            flextable::save_as_docx(
                path = str_glue("{od}/Table_Clinical_Multi_Cox_{x}.docx")
                # pr_section = sect_properties
            )
    }
}

ClinicalCox_DF <- function(clinical = NULL, features_col = c("Smoking", "EGFR", "ALK_eml4", "KRAS", "Histological.Type", "Score"),
                      SaveFile = T, od = "./",p_cut = 2,var_name = NULL) {

    require(tidyverse)
    require(survival)
    require(survminer)
    require(cowplot)
    require(ggpubr)
    require(flextable)
    require(officer)
    require(gtsummary)


    if(is.null(var_name)){
        var_name <- Sys.time() %>% as.character()
    }

    levels <- map(features_col, \(x){
        clinical %>%
            select(x) %>%
            pull() %>%
            factor() %>%
            as.list() %>%
            lvls_union()
    })
    names(levels) <- features_col

    features_chr <- map(features_col, \(x){
        # x <- 'TCGA_Subtype'
        str_c(x, " ", paste0("(", levels[[x]][-1], " vs. ", levels[[x]][1], ")"))
    }) %>%
        unlist() %>%
        as_tibble()


    vars_for_table <- clinical %>%
        dplyr::select(any_of(features_col)) %>%
        colnames(.)
    nv_var <- vars_for_table

    univ_formulas <- sapply(vars_for_table, function(x) as.formula(paste("Surv(time, status)~", x)))
    univ_models <- lapply(univ_formulas, function(x) coxph(x, data = clinical))

    univcox_clinical_res <- map_dfr(.x = univ_models, .f = function(x) {
        result <- summary(x)

        result2 <- data.frame(
            var = rownames(result$coefficients),
            coef = result$coefficients[, 1],
            Hazard_Ratio = result$coefficients[, 2],
            se = result$coefficients[, 3],
            # p.value = result$coefficients[, 5],
            p.value = result$sctest["pvalue"],
            lower_.95 = result$conf.int[, "lower .95"],
            upper_.95 = result$conf.int[, "upper .95"],
            logrank_pvalue = result$sctest["pvalue"],
            wald_pvalue = result$waldtest["pvalue"],
            Likelihood_pvalue = result$logtest["pvalue"]
        )

        result2 <- result2 %>%
            mutate(HR = paste0(
                sprintf("%.2f", Hazard_Ratio),
                "(",
                sprintf("%.2f", lower_.95),
                "-",
                sprintf("%.2f", upper_.95),
                ")"
            ))
    })

    rownames(univcox_clinical_res) <- univcox_clinical_res$var
    univcox_clinical_res$var <- features_chr$value

    univ_cox_df <- univcox_clinical_res %>%
        mutate(P = ifelse(p.value < 0.001, "<0.001", sprintf("%.3f", p.value))) %>%
        select(var, HR, P) %>%
        rename("HR(95%CI)" = HR)

    # MultiCox Part-------------

    # p_cut <- .05
    do_feature_select <- function(p_cut) {
        col_names_left <- map(univ_cox_df %>% filter(P < p_cut) %>% pull(1) %>% str_split(" \\("), ~ .x[[1]]) %>%
            unlist() %>%
            unique()

        vars_for_table <- clinical %>%
            dplyr::select(any_of(col_names_left)) %>%
            colnames(.)

        levels <- map(col_names_left, \(x){
            clinical %>%
                select(x) %>%
                pull() %>%
                factor() %>%
                as.list() %>%
                lvls_union()
        })
        names(levels) <- col_names_left

        features_chr <- map(col_names_left, \(x){
            # x <- 'TCGA_Subtype'
            str_c(x, " ", paste0("(", levels[[x]][-1], " vs. ", levels[[x]][1], ")"))
        }) %>%
            unlist() %>%
            as_tibble()

        return(list("features_chr" = features_chr, "col_names_left" = col_names_left))
    }

    features_new <- do_feature_select(p_cut)

    if (nrow(features_new$features_chr) < 2) {
        features_chr <- features_chr
        vars_for_table <- vars_for_table

        warning(str_glue("在p小于 {crayon::blue(p_cut)} 的时候，只有一个临床信息保留，所以使用全部变量进行多因素Cox分析。"))
    } else {
        features_chr <- features_new$features_chr
        vars_for_table <- features_new$col_names_left
    }


    multicox_formulas <-
        as.formula(paste(
            "Surv(time, status)~",
            paste0(sep = "`", vars_for_table, sep = "`", collapse = "+")
        ))


    # new_style part-----------
    tbl_uv <- tbl_uvregression(
        clinical[c(nv_var, "time", "status")],
        method = coxph,
        y = Surv(time, status),
        exponentiate = T,
        pvalue_fun = function(x) ifelse(x < 0.001, "<0.001", sprintf("%.3f", x))
    ) %>% gtsummary::add_overall()

    tbl_mv <- coxph(
        formula = multicox_formulas,
        data = clinical
    ) %>%
        tbl_regression(
            exponentiate = TRUE, intercept = T,
            pvalue_fun = function(x) ifelse(x < 0.001, "<0.001", sprintf("%.3f", x))
        ) 

    tbl_merge_uv_mv <- tbl_merge(
        tbls = list(tbl_uv, tbl_mv),
        tab_spanner = c("**Univariate Cox**", "**Multivariate Cox**")
    )

    # 设置word格式
    # sect_properties <- prop_section(
    #     page_size = page_size(
    #         orient = "portrait",
    #         width = 21 / 2.54, height = 29.7 / 2.54
    #     ),
    #     type = "continuous",
    #     page_margins = page_mar()
    # )

    if (isTRUE(SaveFile)) {
        if(!dir.exists(od)){
            dir.create(od,recursive = T)
        }

        tbl_merge_uv_mv %>%
            as_hux_xlsx(file = str_glue("{od}/Table_Clinical_Univ_Multi_Cox_{var_name}.xlsx"))
        
        tbl_uv %>%
            as_hux_xlsx(file = str_glue("{od}/Table_Clinical_Univ_Cox_{var_name}.xlsx"))

        tbl_mv %>%
            as_hux_xlsx(file = str_glue("{od}/Table_Clinical_Multi_Cox_{var_name}.xlsx"))        


        tbl_merge_uv_mv %>%
            as_flex_table() %>%
            flextable::save_as_docx(
                path = str_glue("{od}/Table_Clinical_Univ_Multi_Cox_{var_name}.docx")
                # pr_section = sect_properties
            )
        tbl_uv %>%
            as_flex_table() %>%
            flextable::save_as_docx(
                path = str_glue("{od}/Table_Clinical_Univ_Cox_{var_name}.docx")
                # pr_section = sect_properties
            )
        tbl_mv %>%
            as_flex_table() %>%
            flextable::save_as_docx(
                path = str_glue("{od}/Table_Clinical_Multi_Cox_{var_name}.docx")
                # pr_section = sect_properties
            )
    }
}





